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We consider the use of low speed preconditioning for 
time dependent problems. These are solved using a dual 
time step approach. We consider the effect of this dual 
time step on the parameter of the low speed precondi- 
tioning. In addition, we compare the use of two sets of 
variables, conservation and primitive variables, to solve 
the system. We show the effect of these choices on both 
the convergence to a steady state and the accuracy of the 
numerical solutions for low Mach number steady state 
and time dependent flows. 


to the preconditioning matrix. In the next section we 
will discuss in more details the options for choosing w. 
Equation (1) is advanced in artificial time by a Runge- 
Kutta (RK) method until a steady state in r is reached. 
We replace the physical time derivative by a backward 
difference formula (BDF). The general formula for BDF 
schemes can be presented as 

dw c t w n+1 - F(w n ,w n ~ 1 ,...) 
dt At 


Introduction 

Methods for preconditioning the low speed Euler and 
Navier-Stokes equations have been available for about 
twenty years. In spite of these years of development and 
analysis there still exists difficulties with the robustness 
of these techniques. This manifests itself in that some of 
the parameters of the preconditioning matrix are problem 
dependent. Because of the stagnation points it is neces- 
sary to prevent the parameters from becoming too small. 
Unfortunately, the cutoff is frequently large and problem 
dependent. In addition the theory of preconditioning is 
mainly based on properties of the inviscid Euler equa- 
tions with artificial viscosity even though most of the 
applications are for viscous flows. 

The preconditioning changes the time dependent be- 
havior of the system and so is only directly useful for 
steady state calculations. To overcome this difficulty a 
dual time step method has been used by many inves- 
tigators. 3-6,20 In this approach the solution at the next 
physical time step is determined as a steady state prob- 
lem to which preconditioning is applicable. We shall 
analyze time dependent effects on the preconditioning 
parameters. We consider 


P 1 


dw 

df 


dw 

Ht 


+ R = 0 


( 1 ) 


where t is the physical time, r is an artificial time, R 
denotes the residual for the steady Navier-Stokes equa- 
tions, w refers to a general set of unknowns and P refers 
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where c t is constant which depends on the choice of BDF 
scheme. Let superscript 0 denote the last artificial time 
step, k the most recent stage of RK, n the last physical 
time step and n + 1 the next physical time step. A typical 
stage of the RK is 


„ k + 1 _ „„0 


= ur — afcArP < R 


c t w n+1 - F{w n , ...) 
At 


where a k are the stage coefficients of the RK scheme. 
In practice, only the inviscid portion of R k is updated at 
each stage. The viscous portion is updated, for a 5 stage 
scheme, only on the odd stages. Because w n+l is not 
known, we replace it by w k+1 , i.e. current stage of RK. 
We reformulate this as 


w 


fc+i 


= w° — ctfcArP < R 


, fc c t w k - F(w n , ...) 


At 


- aifcCt A tP( 


w k + 1 - w k . 


At 


We apply residual smoothing to the term inside the curly 
brackets. This is done so that the residual smoothing op- 
erates on a difference that vanishes in the steady state. 
Collecting terms we have 


(I + a k c t ^P)w k+1 - ^ 


-a fc ArP j 

At k 
+ ot k c t —Vw 


)W = w 
,k 


c t w k — F(w n , ...) 

A t 


( 2 ) 


The space discretization consists of a central differ- 
ence formula plus a matrix valued artificial dissipation 
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using second and fourth differences. We now describe 
the artificial viscosity due to the second order differenc- 
ing. We express the dissipation in terms of derivatives 
rather than differences for presentation only. We con- 
sider the equation 


dw dF dG 
dt + dx dy 


dH 

— — = vise terms 
dz 


^qAQ, we get 


(I + a k c t ^P Q )Q k+1 = Q° 


- a fc Arr <j R k Q + 


l 


c t w k - F(w ™, ...) 
At 


At . 


+ a kCt~^PQQ 


(6) 


where F, G, and Ft represent the fluxes along x, y, and z 
directions. In the x direction, we supplement the residual 
by the artificial viscosity as follows: 

r =S- , ‘4<« p “ 1 ' pa 'S> ® 

where A is the Jacobian of F with respect to w, h x is the 
mesh spacing, and e 2 is a constant. The absolute value of 
the matrix is found by diagonalizing the matrix, taking 
absolute values of the eigenvalues, with appropriate cut- 
offs to avoid singularities. A similar procedure is used 
for the y and z directions. 

Choice of Variables 

We consider the sets of variables defined by 


w c = (p, pu , pv, pw, E) 

Q = (p,u,v,w,T) (4) 


Wo = ( p,u,v,w,S ) , 


dwo = 


( — , du, dv, dw, dS) 
pc 


We shall refer to Q as the primitive variables and w c as 
the conservation variables. The flux F and the physical 
time derivatives are evaluated in conservation variables 
so that the correct shock jumps are obtained. One im- 
plementation of (2) is to use the conservation variables 
throughout the equation. As the Mach number decreases 
to zero, the density usually becomes constant and so the 
conservation variables become less accurate. For almost 
incompressible flow the primitive variables are more ap- 
propriate. 

If we change from conservation variables to Q vari- 
ables we replace (3) by 


Rq — 
Pq = 
r -1 = 


- h 


dF 
dx 
dQ 
dw c 

dw c-n-i 

dQ Q 


d 

dx 

dw c 

dQ 


(e 2 T~ 1 \P Q A Q \^) 


dx ' 


r = p f 


dQ _ dQ 
dw r dw r 


(5) 


Multiplying (2) by ^ and substituting Aw c = 


After each stage w k+1 is calculated using the nonlinear 
relation between w c and Q. Note that Q only appears in 
the artificial time terms and the artificial viscosity. After 
the artificial time derivative approaches zero the resultant 
equation is in conservation form including the physical 
time derivative. If the physical time derivative would 
also be transformed to Q variables then we might lose 
the conservation form and hence the correct jump condi- 
tions at a shock. Preconditioning destroys conservation 
in the midst of the pseudo-time iteration process. How- 
ever, when the pseudo-time derivative approaches zero, 
the algorithm should recover the conservation form. 

Another possibility is to consider a mixture of conser- 
vation and Q variables. When evaluating the artificial 
viscosity we use Q variables as given by (5). However, 
when updating the variables we revert to w c variables. 
This would be equivalent to (6) if the relation between 
w c and Q variables were linear. We then get 


(/ + a k c t 


At 
A t 


P)w k c +1 


= w 


0 

c 


-a k ArP^R k Q 
At k 

+ CXkCt ~At Pw * 


c t w k 



(7) 


This is the same as (2) except that the artificial viscosity 
in w c variables is replaced by Rq based on Q variables. 
Hence, the two steady states (within the time dependent 
problem) are different while (6) and (7) have the same 
numerical steady state. 

Computations demonstrate that the variables used in 
the artificial viscosity have a much larger effect than the 
choice of variables used to update the solution. Hence, 
we shall concentrate on comparing (2) with (6) and less 
on the mixed formulation (7). We shall further see that 
we can efficiently solve these linear systems. 


Low Speed Preconditioning 

In the above description P is a preconditioning oper- 
ator based on the conservation variables, which is a full 
matrix and is difficult to analyze. Instead, we consider 
the entropy variables vjq ■ Now, the entropy equation 
decouples from the other variables. Furthermore, the Ja- 
cobian matrix is sparse. The simplest preconditioner in 
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w o variables is given by, see. 

9, 10, 12, 14 


where c„ = . Let 
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, (7-1 )T 
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(3 is a parameter which should be of the order of the Mach 
number as to approximately equalize all the eigenvalues 
of P 0 A 0 . 

Let c 2 = ' m , q 2 = u 2 + v 2 + w 2 and q 2 = G-Vq . 

Then the Jacobians that connect these variables are 


then for (6) 


(/ + d-P Q ) 1 x = 


x + ecciC 
1 + d 


d{ 1-/3 2 ) 
1 + d-p 2 


( r 

dwo _ 
dw c 

1-2 P 2 
\q -c 


(1-7 )« 

l 

p 

0 

0 

(l-7)w 


(1 — 7)v (1 — 7)10 7— 1\ 


(l — 7)ti (1 — 7)10 7— 1/ 


dw c 

dw 0 




3* 

h 

, 


pu pv 


0 

0 

0 

p 

pw 


\ 

U 

~ ~71 


'^2 
M 2 / 

2 / 


where h = r 

7-1 


The preconditioner P c in conser- 


dw c 


dw 0 


vation variables is then given by P c = 
calculate P c times a vector x we do it in stages. 


7— 1 

yi = J 72- 


^-27 - (ux 2 + ra 3 + 10x4) + £5 

/A 


and 


z = 


u 

v 

w 

W 


To 


Then 


P c f = x + (/3 2 - 1)2/1 z (8) 

P" 1 ^ = x + - 1)2/1 £ 


In (7) we need to evaluate (/ + d ■ P c ) 1 times a vector 
where d = otkCt^- Then 


(I + d- P c ) 1 x = 


x + eyiZ 
1 + d 


(1 ~P 2 )d 
1 + (3 2 d 


For the primitive variables, Q, we have 


J. 

( ( 3 2 0 

0 

0 

°\ 


0 1 

0 

0 

0 

<0 

II 

0 0 

1 

0 

0 

0 0 

0 

1 

0 

1 

(/3 2 -l )t q 

\ C p p 

0 

0 
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Note,r^ = P Q . 

We need to choose j3 2 and the pseudo-time step. When 
we ignore the correction term and use an explicit formula 
it requires that the pseudo-time step also include a phys- 
ical time step contribution. The amplification matrix in 
pseudo-time for u>o variables, for the two dimensional 
Euler equation in generalized coordinates is given by 

G(ff) = Po — b oj\A + W2-B^ (9) 

where Vol is the volume of the cell and A,B are the Ja- 
cobian matrices of the inviscid flux vectors in the two 
generalized coordinate space dimensions. A and B are 
symmetric in u>o variables and so this is a symmetric hy- 
perbolic system where the physical t derivative is treated 
as another space derivative relative to the marching di- 
rection r. We denote the surface area of the cell as S t) 
where the first subscript refers to the direction of the nor- 
mal and the second is the projection of that normal in 
each direction. Define the contravariant velocity compo- 
nents as 


U uS X x vS xy “b ixS xz 

V = USy X + VSyy + WSy Z ( 10 ) 

W = uS zx + vS zy + wS zz 

In two space dimensions 


G[9) = uj 0 


c t PpAf 

Vol 


u 1 


■ U>2 


/3 2 U 

(3 2 cS xx 

(3 2 cS X y 

C$xx 

U 

0 

C$xy 

0 

u 

'(3 2 V 

(3 2 cS yx 

P 2 cSyy 

cSyx 

V 

0 

r Q 

GUyy 

0 

V 


In three dimensions let 


D = 



c t Vol 

At 


2 

+ U 2 +V 2 +W 2 . 
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Then choose 

c t Vol 
“ 0= ^tD 


u 1 


u 

V 

w 

D 

id 2 — 

D 

^ - D 


and so uiQ+uif +w|+w| = 1. This makes G{8) in (9) a 
convex combination of the Jacobian matrices. Define 

fr S XX U + S yx V + S ZX W 
D- Vol 2 

> _ S X yU + SyyV + S Z yW 
D ■ Vol 2 

s xz u + s yz v + s zz w 

D ■ Vol 2 
q 2 = U 2 + V 2 + W 2 


Thus U, V . W are combinations of the velocity compo- 
nents u, v, w that depend on the geometry metrics. This 
combination comes from the analysis and does not nec- 
essarily have any physical interpretation. This gives 

//3 2 D (3 2 cU f3 2 cV (3 2 cW\ 
cU D 0 0 

( ' cV 0 D 0 
\cW 0 0 D ) 

We can easily symmetrize G but this is not necessary in 
our case. The eigenvalues are A o = D and 


Similarly, when A / phvs is small enough then precondi- 
tioning is turned off globally. 

From (11) we define A i nv = A+. We scale A V i S by 
the volume so it also has dimensions J— . We choose the 

time 

total time step by 


At = 


1 


^inv + uAvis + 


bc t 

At 


(13) 


where a and 6 are constants. Typical values are a = 4—5 
and b = 9 — 25. These are determined semi-empirically 
based on numerical experiments. Straightforward linear 
analysis shows that using a simple implicit formula in 
each stage of the Runge-Kutta formula does not allow 
the use of a larger time step. Hence, even in this case we 
choose b = 25. In order to reduce this value of b one 
would need an explicit-implicit type formula as derived 
in 1 . Note that even when 6 = 0, Ar depends on Af p h ys 
through Ai nv which is a function of (3 2 . However, for the 
explicit treatment of the time derivative we need to add 
a term to the time step that depends on ^ even when 
preconditioning is not used (i.e. (3 = 1). Our choice 
for 6 and Ar is different than that of Venkateswaran 
and Merkle. 18 In addition the parameters of the resid- 
ual smoothing should now depend on the physical time 
derivative term. 


X± = l3 ^D±^(^ J ^j 2 D 2 +/3 2 c 2 q 2 (11) 

The analysis so far is based on inviscid equations: there- 
fore f3 should be designated as /3j ra „, which is chosen 
such that (3cq = Ao = D. For low Mach numbers /3 2 
is small and so A± ~ 1+ 2 V ^ ZJ. We then choose j3 im as a 
term that depends on D plus a cutoff to prevent / 3 from 
becoming too small. This cutoff depends on a global 
quantity M re f. We choose 

0L = K\ + K 2 M 2 f ( 12 ) 

c z q J 

We stress that (3 2 depends on c 0' 1 even when using an 
implicit method for the time derivative term. K\ and 
K 2 are user defined constants. In most cases /\' 2 can 
be chosen as zero since the physical time derivative term 
prevents (3 from becoming too small. This is in contrast 
to steady state problems where this term is crucial. 

We also account for the viscous effects in (3 2 but this 
is only used to reduce the time step rather than directly 
changing the (3 used in the preconditioning or the arti- 
ficial viscosity. 15 For exterior problems we find that in 
the farfield where there are large volumes, the precon- 
ditioning is turned off locally due to the time derivative. 


Results 

The governing equations are solved using a finite vol- 
ume central difference code augmented by a matrix arti- 
ficial viscosity. 2,8 The equations are advanced in pseudo- 
time by a five stage Runge-Kutta scheme accelerated by 
residual smoothing and multigrid. 2, 16 The parameters for 
all the cases are identical with the exception that without 
preconditioning the explicit CFL in the residual smooth- 
ing is 3.75, while with preconditioning CFL-explicit is 
3.25. For all preconditioning cases /3^ in = M0. The 
dual time-stepping uses a second order BDF formula 
stabilized with the algorithm of Melson and Sanetrik. 6 
Early work on dual stepping and preconditioning is pre- 
sented in. 7,20 Turkel et al. 13 compared results from two 
different preconditioning schemes for steady state prob- 
lems. 

NACA4412 

We first consider a steady state case, turbulent flow 
around a NACA4412 airfoil. A mesh with 257 x 81 
grid points, constructed using Wigton’s 19 grid genera- 
tion procedure and displayed in Fig. 1, is used for these 
computations. The inflow conditions are a = 13.87° 
and Moo = 0.2, 0.05 and 0.01. The Reynolds number is 
1.52 x 10 6 and the Baldwin-Lomax turbulence model is 
used. 
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We compare the convergence rate without precon- 
ditioning versus preconditioning based on conservation 
variables and primitive (p,u,v,w,T) variables. For all 
preconditioned cases we use K 2 = 1.0. We show a por- 
tion of the grid in figure 1. Every second grid line is 
shown to increase the clarity. There are 50 iterations on 
two coarser grids and 500 five stage RK iterations on the 
finest grid. In figures 2 and 3 we compare the residual as 
well as the drag history for = 0.2. We stress that the 
only difference between the two preconditioned results 
is the artificial viscosity, (Eqns. 2 and 7). In one case 
it is evaluated directly on the conservation variables. In 
the other case we take differences of the primitive vari- 
ables, multiply by r _ 1 |PH| in primitive variables and 
transform back to conservation variables. We see in the 
following figures that in the beginning of the computa- 
tion both approaches give the same rate of convergence 
and accuracy. However, for small residuals an artifi- 
cial viscosity based on the conservation variables stalls. 
When the artificial viscosity is based on the primitive 
variables convergence continues to decrease further. Our 
explanation for this phenomena is that for the conserva- 
tion variables the density is fairly constant. Hence, at low 
residual values the contribution from density derivatives 
is negligible. For the first few iterations the precondition- 
ing slows the convergence rate of the residual. However, 
the drag converges faster even at the beginning using 
preconditioning. Overall convergence is improved sig- 
nificantly when preconditioning is used. 



Fig. 3 NACA4412 - M ^ = 0.2, drag history 


We next consider the same case but with an inflow 
Mach number 0.05. The results are given in figures 4 
- 6. The improved convergence for the precondition- 
ing is more evident for this case. The convergence for 
the conservation variables is reasonable until it bottoms 
out. The convergence with the primitive variables con- 
tinues in a straight line. Note the difference in the fi- 
nal lift coefficient between the preconditioned and non- 
preconditioned schemes. This is a reflection of the de- 
crease in accuracy without preconditioning for low Mach 
numbers. 11, 14 
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Fig. 4 NACA4412 - = 0.05, residual history 



Fig. 5 NACA4412 - M = 0.05, drag history 



CYCLES 

Fig. 6 NACA4412 - M x = 0.05, lift history 


We finally show the results in figures 7 through 9 for 
the same case but an inflow Mach number 0.01. The dif- 
ference in lift and drag between the preconditioned and 
non-preconditioned algorithms is now quite noticeable. 
The pattern of the convergence history is similar to the 
previous cases but is more dramatic. We now clearly see 
that using conservation variables in the artificial viscosity 
limits the convergence of the residual even with precon- 
ditioning. 




0 200 400 600 800 1000 

CYCLES 

Fig. 8 NACA4412 - = 0.01, drag history 
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Fig. 9 NACA4412 - M x = 0.01, lift history 


NACA0012- high angle of attack 


We next consider turbulent flow around a NACA0012 
airfoil. An O mesh containing 141 x 61 points (Fig. 10) is 
used for this configuration. An exploded view of the grid 
near the trailing edge is shown in Fig. 11. The inflow 
conditions are M 0 0 = 0.1, an angle of attack of 12° and 
Re = 3 x 10 6 . A Spalart-Allmaras turbulence model is 
used. At this angle of attack the flow is still steady. We 
also consider an angle of attack of 30° where the flow is 
no longer steady. 



Fig. 10 Partial view of grid for NACA0012 


Fig. 11 Grid near trailing edge of NACA0012 

We begin with the steady state calculation for a = 
12°. The convergence of the residual and drag coefficient 
are shown in Figs. 12 and 13. We see that without pre- 
conditioning the convergence stalls, as expected, with an 
inflow Mach number = 0.1. Basing the precondi- 
tioning on the conservation variables (Eqn. 2) improves 
the convergence rate. However, using primitive variables 
in the artificial viscosity (Eqn. 6) yields a better asymp- 
totic rate. Although not shown here, we observed that 
the variable used for the update ( Eqn. 6 versus Eqn. 
7) has no effect on the convergence rate. When using 
w c variables the residual is based on the density while 
when using Q variables the residual is based on the pres- 
sure. We see that the preconditioning yields significantly 
better convergence than the non-preconditioned scheme. 
Furthermore, the steady state is different when precondi- 
tioning is used. Based on earlier work of Turkel et ah, 11 
we expect the preconditioned results to be more accurate. 

We next examine the convergence and accuracy for the 
time dependent case. We consider an angle of attack with 
a = 30° and inflow Mach number of = 0.1. Since 
the flow is now time dependent we are only interested in 
the convergence rate within a physical time cycle. We 
use 30 cycles (5 stage Runge-Kutta with multigrid and 
residual smoothing) of pseudo-time stepping within each 
physical time cycle. In Fig. 14 we display the residual 
(density for w c variables and pressure for Q variables) 
while in Fig. 15 we display the lift. In both cases we see 
that initially the convergence is best without any precon- 
ditioning. The residual stalls after a few cycles but the 
lift has already reached its new level to within graphical 
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accuracy. The overall convergence rate for the residual 
is improved with preconditioning on conservative and 
primitive variables. The convergence of lift with precon- 
ditioning with an artificial viscosity based on conserva- 
tion, w c variables is slightly worse but is also reasonable. 
When the preconditioning uses primitive variables in the 
artificial viscosity then the initial convergence is slowed 
and one requires about 15-20 subiteration cycles for the 
lift to converge to graphical accuracy. In all cases the 
variables used to advance the solution to the next cycle 
had little effect (figure not shown here) on either the con- 
vergence rate of the subiterations or the on the accuracy. 
The major influence is due to the variables used in the 
artificial viscosity. 


We stress that the use of preconditioning affects the 
accuracy of the solution and so changes the values of the 
lift and drag. In Fig. 16 we plot the time dependent his- 
tory of the lift. We clearly see that for short times the two 
preconditioned results, with the artificial viscosity based 
on conservation or primitive variables, agree with each 
other but give a different lift than the non-preconditioned 
algorithm. Based on the results from low speed steady 
flows that demonstrated improved accuracy with precon- 
ditioning, 11 it is reasonable to assume that the unsteady 
solutions obtained here with preconditioning are more 
accurate compared to the un-preconditioned solutions. 




Fig. 13 Cd for steady state NACA0012 



CYCLES 


Fig. 14 Residual within subiterations for a = 30° 



CYCLES 

Fig. 15 c L within subiterations for a = 30° 
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Fig. 16 Time dependent cl for a = 30° 


Concluding Remarks and Future Directions 

A general preconditioning formulation for treating low 
speed flows with compressible flow equations is pre- 
sented. Several alternative forms of preconditioning vari- 
ables have been examined. The choice of the variables 
used to advance the solution in time had very little effect 
on convergence or accuracy. The convergence for low 
speed steady state problems is improved significantly 
when preconditioning is used. The primitive variables 
based preconditioner appears to be most effective for 
steady state flows. However, for time-dependent prob- 
lems, preliminary results about the choice of precondi- 
tioners are somewhat inconclusive, and requires further 
testing. 
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